/* Copyright (c) 2021, National University of Defense Technology. All rights reserved.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 * 
 *     http://www.apache.org/licenses/LICENSE-2.0
 * 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <set>
#include <map>
#include <stdlib.h>
#include "parcsr_ilu_preconditioner.h"

namespace YHAMG
{

ParCSR_ILU_Preconditioner::ParCSR_ILU_Preconditioner()
	: interior(0),
	label(0),
	num_neighbors(0),
	neighbor_labels(0),
	neighbor_ranks(0),
	recv_ptr(0),
	send_ptr(0),
	send_list(0),
	send_buffer(0)
{
}

ParCSR_ILU_Preconditioner::~ParCSR_ILU_Preconditioner()
{
	if (interior) delete[] interior;
	if (neighbor_labels) delete[] neighbor_labels;
	if (recv_ptr) delete[] recv_ptr;
	if (send_ptr) delete[] send_ptr;
	if (send_list) delete[] send_list;
	if (send_buffer) delete[] send_buffer;
}

void ParCSR_ILU_Preconditioner::Free()
{
	D_Recip_Loc.Free();
	L_LocLoc.Free();
	L_LocExt.Free();
	U_LocLoc.Free();
	U_LocExt.Free();
	x_Ext.Free();
	if (interior) delete[] interior;
	if (neighbor_labels) delete[] neighbor_labels;
	if (recv_ptr) delete[] recv_ptr;
	if (send_ptr) delete[] send_ptr;
	if (send_list) delete[] send_list;
	if (send_buffer) delete[] send_buffer;
	label = 0;
	num_neighbors = 0;
	interior = 0;
	neighbor_labels = 0;
	recv_ptr = 0;
	send_ptr = 0;
	send_list = 0;
	send_buffer = 0;
}

#define MPI_TAG 500

static inline void Swap(int& a, int& b)
{
	int temp = a;
	a = b;
	b = temp;
}

#define CMP(a, b) ((a) < (b))
#define CMPLOC(i, j) (interior[i] != interior[j] ? interior[i] : (i) < (j))
#define CMPEXT(i, j) (from[i] == from[j] ? map[i] < map[j] : label[from[i]] < label[from[j]])

static void Sort(int* a, int left, int right)
{
	if (left >= right) return;

	Swap(a[left], a[(left + right) / 2]);

	int last = left;
	for (int i = left + 1; i <= right; ++i)
		if (CMP(a[i], a[left]))
			Swap(a[++last], a[i]);

	Swap(a[left], a[last]);

	Sort(a, left, last - 1);
	Sort(a, last + 1, right);
}

static void ILUSortLoc(int* a, int left, int right, const bool* interior)
{
	if (left >= right) return;

	Swap(a[left], a[(left + right) / 2]);

	int last = left;
	for (int i = left + 1; i <= right; ++i)
		if (CMPLOC(a[i], a[left]))
			Swap(a[++last], a[i]);

	Swap(a[left], a[last]);

	ILUSortLoc(a, left, last - 1, interior);
	ILUSortLoc(a, last + 1, right, interior);
}

static void ILUSortExt(int* a, int left, int right, const int* map, const int* from, const int* label)
{
	if (left >= right) return;

	Swap(a[left], a[(left + right) / 2]);

	int last = left;
	for (int i = left + 1; i <= right; ++i)
		if (CMPEXT(a[i], a[left]))
			Swap(a[++last], a[i]);

	Swap(a[left], a[last]);

	ILUSortExt(a, left, last - 1, map, from, label);
	ILUSortExt(a, last + 1, right, map, from, label);
}

void ParCSR_ILU_Preconditioner::Setup(const ParCSR_Matrix& A, int REUSE)
{
	int n = A.LocLoc.size[0];
	int Ext_size = A.LocExt.size[1];

	int* A_LocLoc_rowptr = A.LocLoc.rowptr;
	int* A_LocLoc_colind = A.LocLoc.colind;
	double* A_LocLoc_values = A.LocLoc.values;
	int* A_LocExt_rowptr = A.LocExt.rowptr;
	int* A_LocExt_colind = A.LocExt.colind;
	double* A_LocExt_values = A.LocExt.values;

	int* Ext_map = A.Ext_map;
		
	if (REUSE == 0)
	{
		Free();
		
		comm = A.comm;

		int rank;
		MPI_Comm_rank(comm, &rank);
		
		num_neighbors = A.num_neighbors;
		neighbor_ranks = new int[num_neighbors];
		recv_ptr = new int[num_neighbors + 1];
		send_ptr = new int[num_neighbors + 1];
		send_list = new int[A.send_ptr[num_neighbors]];
		send_buffer = new double[A.send_ptr[num_neighbors]];

		for (int r = 0; r < num_neighbors; ++r)
			neighbor_ranks[r] = A.neighbor_ranks[r];
		for (int r = 0; r <= num_neighbors; ++r)
			recv_ptr[r] = A.Ext_ptr[r];
		for (int r = 0; r <= num_neighbors; ++r)
			send_ptr[r] = A.send_ptr[r];
		for (int i = 0; i < A.send_ptr[num_neighbors]; ++i)
			send_list[i] = A.send_list[i];

		int* w = new int[n];
		double* x = new double[n];

		for (int i = 0; i < n; ++i)
			w[i] = -1;

		interior = new bool[n];
		for (int i = 0; i < n; ++i)
			interior[i] = A_LocExt_rowptr[i + 1] == A_LocExt_rowptr[i];
		for (int i = 0; i < send_ptr[num_neighbors]; ++i)
			interior[send_list[i]] = 0;

		MPI_Request* recv_requests = new MPI_Request[num_neighbors];
		MPI_Request* send_requests = new MPI_Request[num_neighbors];
		MPI_Status status;

		label = 0;
		neighbor_labels = new int[num_neighbors];

		std::set<int> set_neighbor_labels;

		for (int r = 0; r < num_neighbors; ++r)
			MPI_Irecv(neighbor_labels + r, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (rank > neighbor_ranks[r])
			{
				MPI_Wait(recv_requests + r, &status);
				set_neighbor_labels.insert(neighbor_labels[r]);
			}
		}

		for (std::set<int>::iterator iter = set_neighbor_labels.begin(); iter != set_neighbor_labels.end(); ++iter, ++label)
			if (*iter != label) break;

		for (int r = 0; r < num_neighbors; ++r)
			MPI_Send(&label, 1, MPI_INT, neighbor_ranks[r], MPI_TAG, comm);

		for (int r = 0; r < num_neighbors; ++r)
			if (rank < neighbor_ranks[r])
				MPI_Wait(recv_requests + r, &status);

		std::map<int, int> map_neighbor;
		std::map<int, int>* map_Ext = new std::map<int, int>[num_neighbors];

		for (int r = 0; r < num_neighbors; ++r)
		{
			map_neighbor.insert(std::pair<int, int>(neighbor_ranks[r], r));
			for (int i = recv_ptr[r]; i < recv_ptr[r + 1]; ++i)
				map_Ext[r].insert(std::pair<int, int>(Ext_map[i], i));
		}

		int* Ext_from = new int[Ext_size];
		for (int r = 0; r < num_neighbors; ++r)
			for (int i = recv_ptr[r]; i < recv_ptr[r + 1]; ++i)
				Ext_from[i] = r;

		D_Recip_Loc.Resize(n);

		double* D_Recip_Loc_values = D_Recip_Loc.values;

		int* L_LocLoc_rowptr = new int[n + 1];
		int* U_LocLoc_rowptr = new int[n + 1];
		int* L_LocExt_rowptr = new int[n + 1];
		int* U_LocExt_rowptr = new int[n + 1];

		L_LocLoc_rowptr[0] = 0;
		U_LocLoc_rowptr[0] = 0;
		L_LocExt_rowptr[0] = 0;
		U_LocExt_rowptr[0] = 0;

		for (int i = 0, cnt_l_loc = 0, cnt_u_loc = 0, cnt_l_ext = 0, cnt_u_ext = 0; i < n; L_LocLoc_rowptr[++i] = cnt_l_loc, U_LocLoc_rowptr[i] = cnt_u_loc, L_LocExt_rowptr[i] = cnt_l_ext, U_LocExt_rowptr[i] = cnt_u_ext)
		{
			if (interior[i])
			{
				for (int j = A_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
				{
					int jcol = A_LocLoc_colind[j];

					if (jcol < i && interior[jcol])
						++cnt_l_loc;
					else if (jcol != i)
						++cnt_u_loc;
				}
			}
			else
			{
				for (int j = A_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
				{
					int jcol = A_LocLoc_colind[j];

					if (jcol < i || interior[jcol])
						++cnt_l_loc;
					else if (jcol != i)
						++cnt_u_loc;
				}

				for (int j = A_LocExt_rowptr[i]; j < A_LocExt_rowptr[i + 1]; ++j)
				{
					int jcol = A_LocExt_colind[j];
					int jfrom = Ext_from[jcol];

					if (label > neighbor_labels[jfrom])
						++cnt_l_ext;
					else
						++cnt_u_ext;
				}
			}
		}

		int* L_LocLoc_colind = new int[L_LocLoc_rowptr[n]];
		double* L_LocLoc_values = new double[L_LocLoc_rowptr[n]];
		int* U_LocLoc_colind = new int[U_LocLoc_rowptr[n]];
		double* U_LocLoc_values = new double[U_LocLoc_rowptr[n]];
		int* L_LocExt_colind = new int[L_LocExt_rowptr[n]];
		double* L_LocExt_values = new double[L_LocExt_rowptr[n]];
		int* U_LocExt_colind = new int[U_LocExt_rowptr[n]];
		double* U_LocExt_values = new double[U_LocExt_rowptr[n]];

		for (int i = 0; i < n; ++i)
		{
			if (interior[i])
			{
				w[i] = i;
				x[i] = 0.0;

				for (int j = A_LocLoc_rowptr[i], k = L_LocLoc_rowptr[i], r = U_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
				{
					int jcol = A_LocLoc_colind[j];
					double jval = A_LocLoc_values[j];

					w[jcol] = i;
					x[jcol] = jval;

					if (jcol < i && interior[jcol])
						L_LocLoc_colind[k++] = jcol;
					else if (jcol != i)
						U_LocLoc_colind[r++] = jcol;
				}

				Sort(L_LocLoc_colind, L_LocLoc_rowptr[i], L_LocLoc_rowptr[i + 1] - 1);

				for (int k = L_LocLoc_rowptr[i]; k < L_LocLoc_rowptr[i + 1]; ++k)
				{
					int kcol = L_LocLoc_colind[k];
					double kval = x[kcol] * D_Recip_Loc_values[kcol];

					x[kcol] = kval;

					for (int j = U_LocLoc_rowptr[kcol]; j < U_LocLoc_rowptr[kcol + 1]; ++j)
					{
						int jcol = U_LocLoc_colind[j];
						double jval = U_LocLoc_values[j];

						if (w[jcol] == i)
							x[jcol] -= kval * jval;
					}
				}

				D_Recip_Loc_values[i] = 1.0 / x[i];
				for (int j = L_LocLoc_rowptr[i]; j < L_LocLoc_rowptr[i + 1]; ++j)
					L_LocLoc_values[j] = x[L_LocLoc_colind[j]];
				for (int j = U_LocLoc_rowptr[i]; j < U_LocLoc_rowptr[i + 1]; ++j)
					U_LocLoc_values[j] = x[U_LocLoc_colind[j]];
			}
		}

		char** recv_buffer = new char* [num_neighbors];
		char** send_buffer = new char* [num_neighbors];

		int Ext_nnz = 0;

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label > neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			{
				int bufsize;
				MPI_Probe(neighbor_ranks[r], MPI_TAG, comm, &status);
				MPI_Get_count(&status, MPI_PACKED, &bufsize);
				recv_buffer[r] = (char*)malloc(bufsize);
				MPI_Irecv(recv_buffer[r], bufsize, MPI_PACKED, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

				int num_rows = recv_ptr[r + 1] - recv_ptr[r];
				int num_nonzeros = (bufsize - num_rows * (sizeof(int) + sizeof(double))) / (2 * sizeof(int) + sizeof(double));
				Ext_nnz += num_nonzeros;
			}
		}

		double* D_Recip_Ext_values = new double[Ext_size];
		int* U_ExtLoc_rowptr = new int[Ext_size + 1];
		int* U_ExtLoc_colind = new int[Ext_nnz];
		double* U_ExtLoc_values = new double[Ext_nnz];
		int* U_ExtExt_rowptr = new int[Ext_size + 1];
		int* U_ExtExt_colind = new int[Ext_nnz];
		double* U_ExtExt_values = new double[Ext_nnz];

		U_ExtLoc_rowptr[0] = 0;
		U_ExtExt_rowptr[0] = 0;

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label > neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			{
				MPI_Wait(recv_requests + r, &status);

				int bufsize;
				MPI_Get_count(&status, MPI_PACKED, &bufsize);

				int position = 0;

				MPI_Unpack(recv_buffer[r], bufsize, &position, D_Recip_Ext_values + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_DOUBLE, comm);

				for (int i = recv_ptr[r], k = U_ExtLoc_rowptr[i], t = U_ExtExt_rowptr[i]; i < recv_ptr[r + 1]; U_ExtLoc_rowptr[++i] = k, U_ExtExt_rowptr[i] = t)
				{
					int length;
					MPI_Unpack(recv_buffer[r], bufsize, &position, &length, 1, MPI_INT, comm);

					int* colfrom = new int[length];
					int* colind = new int[length];
					double* values = new double[length];

					MPI_Unpack(recv_buffer[r], bufsize, &position, colfrom, length, MPI_INT, comm);
					MPI_Unpack(recv_buffer[r], bufsize, &position, colind, length, MPI_INT, comm);
					MPI_Unpack(recv_buffer[r], bufsize, &position, values, length, MPI_DOUBLE, comm);

					for (int j = 0; j < length; ++j)
					{
						if (colfrom[j] == rank)
						{
							U_ExtLoc_colind[k] = colind[j];
							U_ExtLoc_values[k++] = values[j];
						}
						else
						{
							std::map<int, int>::iterator ret1 = map_neighbor.find(colfrom[j]);
							if (ret1 == map_neighbor.end()) continue;
							std::map<int, int>::iterator ret2 = map_Ext[ret1->second].find(colind[j]);
							if (ret2 == map_Ext[ret1->second].end()) continue;
							U_ExtExt_colind[t] = ret2->second;
							U_ExtExt_values[t++] = values[j];
						}
					}

					delete[] colfrom;
					delete[] colind;
					delete[] values;
				}

				free(recv_buffer[r]);
			}
			else if (recv_ptr[r + 1] > recv_ptr[r])
			{
				for (int i = recv_ptr[r], k = U_ExtLoc_rowptr[i], t = U_ExtExt_rowptr[i]; i < recv_ptr[r + 1]; U_ExtLoc_rowptr[++i] = k, U_ExtExt_rowptr[i] = t)
					;
			}
		}

		delete[] map_Ext;

		int* v = new int[Ext_size];
		double* y = new double[Ext_size];

		for (int i = 0; i < Ext_size; ++i)
			v[i] = -1;

		for (int i = 0; i < n; ++i)
		{
			if (interior[i]) continue;

			w[i] = i;
			x[i] = 0.0;

			for (int j = A_LocLoc_rowptr[i], k = L_LocLoc_rowptr[i], r = U_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
			{
				int jcol = A_LocLoc_colind[j];
				double jval = A_LocLoc_values[j];

				w[jcol] = i;
				x[jcol] = jval;

				if (jcol < i || interior[jcol])
					L_LocLoc_colind[k++] = jcol;
				else if (jcol != i)
					U_LocLoc_colind[r++] = jcol;
			}

			ILUSortLoc(L_LocLoc_colind, L_LocLoc_rowptr[i], L_LocLoc_rowptr[i + 1] - 1, interior);

			for (int j = A_LocExt_rowptr[i], k = L_LocExt_rowptr[i], r = U_LocExt_rowptr[i]; j < A_LocExt_rowptr[i + 1]; ++j)
			{
				int jcol = A_LocExt_colind[j];
				double jval = A_LocExt_values[j];
				int jfrom = Ext_from[jcol];

				v[jcol] = i;
				y[jcol] = jval;

				if (label > neighbor_labels[jfrom])
					L_LocExt_colind[k++] = jcol;
				else
					U_LocExt_colind[r++] = jcol;
			}

			ILUSortExt(L_LocExt_colind, L_LocExt_rowptr[i], L_LocExt_rowptr[i + 1] - 1, Ext_map, Ext_from, neighbor_labels);

			for (int k = L_LocExt_rowptr[i]; k < L_LocExt_rowptr[i + 1]; ++k)
			{
				int kcol = L_LocExt_colind[k];
				double kval = y[kcol] * D_Recip_Ext_values[kcol];

				y[kcol] = kval;

				for (int j = U_ExtExt_rowptr[kcol]; j < U_ExtExt_rowptr[kcol + 1]; ++j)
				{
					int jcol = U_ExtExt_colind[j];
					double jval = U_ExtExt_values[j];

					if (v[jcol] == i)
						y[jcol] -= kval * jval;
				}

				for (int j = U_ExtLoc_rowptr[kcol]; j < U_ExtLoc_rowptr[kcol + 1]; ++j)
				{
					int jcol = U_ExtLoc_colind[j];
					double jval = U_ExtLoc_values[j];

					if (w[jcol] == i)
						x[jcol] -= kval * jval;
				}
			}

			for (int k = L_LocLoc_rowptr[i]; k < L_LocLoc_rowptr[i + 1]; ++k)
			{
				int kcol = L_LocLoc_colind[k];
				double kval = x[kcol] * D_Recip_Loc_values[kcol];

				x[kcol] = kval;

				for (int j = U_LocLoc_rowptr[kcol]; j < U_LocLoc_rowptr[kcol + 1]; ++j)
				{
					int jcol = U_LocLoc_colind[j];
					double jval = U_LocLoc_values[j];

					if (w[jcol] == i)
						x[jcol] -= kval * jval;
				}
			}

			D_Recip_Loc_values[i] = 1.0 / x[i];
			for (int j = L_LocLoc_rowptr[i]; j < L_LocLoc_rowptr[i + 1]; ++j)
				L_LocLoc_values[j] = x[L_LocLoc_colind[j]];
			for (int j = U_LocLoc_rowptr[i]; j < U_LocLoc_rowptr[i + 1]; ++j)
				U_LocLoc_values[j] = x[U_LocLoc_colind[j]];
			for (int j = L_LocExt_rowptr[i]; j < L_LocExt_rowptr[i + 1]; ++j)
				L_LocExt_values[j] = y[L_LocExt_colind[j]];
			for (int j = U_LocExt_rowptr[i]; j < U_LocExt_rowptr[i + 1]; ++j)
				U_LocExt_values[j] = y[U_LocExt_colind[j]];
		}

		delete[] w;
		delete[] x;
		delete[] v;
		delete[] y;

		delete[] D_Recip_Ext_values;
		delete[] U_ExtLoc_rowptr;
		delete[] U_ExtLoc_colind;
		delete[] U_ExtLoc_values;
		delete[] U_ExtExt_rowptr;
		delete[] U_ExtExt_colind;
		delete[] U_ExtExt_values;

		double* D_Recip_Send = A.send_buffer;

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label < neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
			{
				int num_rows = send_ptr[r + 1] - send_ptr[r];
				int num_nonzeros = 0;
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					num_nonzeros += U_LocLoc_rowptr[send_list[i] + 1] - U_LocLoc_rowptr[send_list[i]] + U_LocExt_rowptr[send_list[i] + 1] - U_LocExt_rowptr[send_list[i]];

				int bufsize = (num_rows + 2 * num_nonzeros) * sizeof(int) + (num_rows + num_nonzeros) * sizeof(double);
				send_buffer[r] = (char*)malloc(bufsize);

				int position = 0;

				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					D_Recip_Send[i] = D_Recip_Loc_values[send_list[i]];

				MPI_Pack(D_Recip_Send + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_DOUBLE, send_buffer[r], bufsize, &position, comm);

				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				{
					int irow = send_list[i];
					int length = U_LocLoc_rowptr[irow + 1] - U_LocLoc_rowptr[irow] + U_LocExt_rowptr[irow + 1] - U_LocExt_rowptr[irow];

					int* colfrom = new int[length];
					int* colind = new int[length];
					double* values = new double[length];

					int j = 0;
					for (int k = U_LocLoc_rowptr[irow]; k < U_LocLoc_rowptr[irow + 1]; ++k)
					{
						colfrom[j] = rank;
						colind[j] = U_LocLoc_colind[k];
						values[j++] = U_LocLoc_values[k];
					}
					for (int k = U_LocExt_rowptr[irow]; k < U_LocExt_rowptr[irow + 1]; ++k)
					{
						int kcol = U_LocExt_colind[k];
						colfrom[j] = neighbor_ranks[Ext_from[kcol]];
						colind[j] = Ext_map[kcol];
						values[j++] = U_LocExt_values[k];
					}

					MPI_Pack(&length, 1, MPI_INT, send_buffer[r], bufsize, &position, comm);
					MPI_Pack(colfrom, length, MPI_INT, send_buffer[r], bufsize, &position, comm);
					MPI_Pack(colind, length, MPI_INT, send_buffer[r], bufsize, &position, comm);
					MPI_Pack(values, length, MPI_DOUBLE, send_buffer[r], bufsize, &position, comm);

					delete[] colfrom;
					delete[] colind;
					delete[] values;
				}

				MPI_Isend(send_buffer[r], bufsize, MPI_PACKED, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label < neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
			{
				MPI_Wait(send_requests + r, &status);
				free(send_buffer[r]);
			}
		}

		delete[] recv_buffer;
		delete[] recv_requests;
		delete[] send_buffer;
		delete[] send_requests;

		delete[] Ext_from;

		L_LocLoc.size[0] = n;
		L_LocLoc.size[1] = n;
		L_LocLoc.rowptr = L_LocLoc_rowptr;
		L_LocLoc.colind = L_LocLoc_colind;
		L_LocLoc.values = L_LocLoc_values;

		L_LocExt.size[0] = n;
		L_LocExt.size[1] = Ext_size;
		L_LocExt.rowptr = L_LocExt_rowptr;
		L_LocExt.colind = L_LocExt_colind;
		L_LocExt.values = L_LocExt_values;

		U_LocLoc.size[0] = n;
		U_LocLoc.size[1] = n;
		U_LocLoc.rowptr = U_LocLoc_rowptr;
		U_LocLoc.colind = U_LocLoc_colind;
		U_LocLoc.values = U_LocLoc_values;

		U_LocExt.size[0] = n;
		U_LocExt.size[1] = Ext_size;
		U_LocExt.rowptr = U_LocExt_rowptr;
		U_LocExt.colind = U_LocExt_colind;
		U_LocExt.values = U_LocExt_values;

		x_Ext.Resize(Ext_size);
	}
	else
	{
		int rank;
		MPI_Comm_rank(comm, &rank);

		int* L_LocLoc_rowptr = L_LocLoc.rowptr;
		int* L_LocLoc_colind = L_LocLoc.colind;
		double* L_LocLoc_values = L_LocLoc.values;
		int* U_LocLoc_rowptr = U_LocLoc.rowptr;
		int* U_LocLoc_colind = U_LocLoc.colind;
		double* U_LocLoc_values = U_LocLoc.values;
		int* L_LocExt_rowptr = L_LocExt.rowptr;
		int* L_LocExt_colind = L_LocExt.colind;
		double* L_LocExt_values = L_LocExt.values;
		int* U_LocExt_rowptr = U_LocExt.rowptr;
		int* U_LocExt_colind = U_LocExt.colind;
		double* U_LocExt_values = U_LocExt.values;
		double* D_Recip_Loc_values = D_Recip_Loc.values;

		int* w = new int[n];
		double* x = new double[n];
		for (int i = 0; i < n; ++i)
			w[i] = -1;

		for (int i = 0; i < n; ++i)
		{
			if (interior[i])
			{
				w[i] = i;
				x[i] = 0.0;

				for (int j = A_LocLoc_rowptr[i], k = L_LocLoc_rowptr[i], r = U_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
				{
					int jcol = A_LocLoc_colind[j];
					double jval = A_LocLoc_values[j];

					w[jcol] = i;
					x[jcol] = jval;
				}

				for (int k = L_LocLoc_rowptr[i]; k < L_LocLoc_rowptr[i + 1]; ++k)
				{
					int kcol = L_LocLoc_colind[k];
					double kval = x[kcol] * D_Recip_Loc_values[kcol];

					x[kcol] = kval;

					for (int j = U_LocLoc_rowptr[kcol]; j < U_LocLoc_rowptr[kcol + 1]; ++j)
					{
						int jcol = U_LocLoc_colind[j];
						double jval = U_LocLoc_values[j];

						if (w[jcol] == i)
							x[jcol] -= kval * jval;
					}
				}

				D_Recip_Loc_values[i] = 1.0 / x[i];
				for (int j = L_LocLoc_rowptr[i]; j < L_LocLoc_rowptr[i + 1]; ++j)
					L_LocLoc_values[j] = x[L_LocLoc_colind[j]];
				for (int j = U_LocLoc_rowptr[i]; j < U_LocLoc_rowptr[i + 1]; ++j)
					U_LocLoc_values[j] = x[U_LocLoc_colind[j]];
			}
		}

		MPI_Request* recv_requests = new MPI_Request[num_neighbors];
		MPI_Request* send_requests = new MPI_Request[num_neighbors];
		MPI_Status status;

		std::map<int, int> map_neighbor;
		std::map<int, int>* map_Ext = new std::map<int, int>[num_neighbors];

		for (int r = 0; r < num_neighbors; ++r)
		{
			map_neighbor.insert(std::pair<int, int>(neighbor_ranks[r], r));
			for (int i = recv_ptr[r]; i < recv_ptr[r + 1]; ++i)
				map_Ext[r].insert(std::pair<int, int>(Ext_map[i], i));
		}

		int* Ext_ranks = new int[Ext_size];
		for (int r = 0; r < num_neighbors; ++r)
			for (int i = recv_ptr[r]; i < recv_ptr[r + 1]; ++i)
				Ext_ranks[i] = neighbor_ranks[r];

		char** recv_buffer = new char* [num_neighbors];
		char** send_buffer = new char* [num_neighbors];

		int Ext_nnz = 0;

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label > neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			{
				int bufsize;
				MPI_Probe(neighbor_ranks[r], MPI_TAG, comm, &status);
				MPI_Get_count(&status, MPI_PACKED, &bufsize);
				recv_buffer[r] = (char*)malloc(bufsize);
				MPI_Irecv(recv_buffer[r], bufsize, MPI_PACKED, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

				int num_rows = recv_ptr[r + 1] - recv_ptr[r];
				int num_nonzeros = (bufsize - num_rows * (sizeof(int) + sizeof(double))) / (2 * sizeof(int) + sizeof(double));
				Ext_nnz += num_nonzeros;
			}
		}

		double* D_Recip_Ext_values = new double[Ext_size];
		int* U_ExtLoc_rowptr = new int[Ext_size + 1];
		int* U_ExtLoc_colind = new int[Ext_nnz];
		double* U_ExtLoc_values = new double[Ext_nnz];
		int* U_ExtExt_rowptr = new int[Ext_size + 1];
		int* U_ExtExt_colind = new int[Ext_nnz];
		double* U_ExtExt_values = new double[Ext_nnz];

		U_ExtLoc_rowptr[0] = 0;
		U_ExtExt_rowptr[0] = 0;

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label > neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			{
				MPI_Wait(recv_requests + r, &status);

				int bufsize;
				MPI_Get_count(&status, MPI_PACKED, &bufsize);

				int position = 0;

				MPI_Unpack(recv_buffer[r], bufsize, &position, D_Recip_Ext_values + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_DOUBLE, comm);

				for (int i = recv_ptr[r], k = U_ExtLoc_rowptr[i], t = U_ExtExt_rowptr[i]; i < recv_ptr[r + 1]; U_ExtLoc_rowptr[++i] = k, U_ExtExt_rowptr[i] = t)
				{
					int length;
					MPI_Unpack(recv_buffer[r], bufsize, &position, &length, 1, MPI_INT, comm);

					int* colfrom = new int[length];
					int* colind = new int[length];
					double* values = new double[length];

					MPI_Unpack(recv_buffer[r], bufsize, &position, colfrom, length, MPI_INT, comm);
					MPI_Unpack(recv_buffer[r], bufsize, &position, colind, length, MPI_INT, comm);
					MPI_Unpack(recv_buffer[r], bufsize, &position, values, length, MPI_DOUBLE, comm);

					for (int j = 0; j < length; ++j)
					{
						if (colfrom[j] == rank)
						{
							U_ExtLoc_colind[k] = colind[j];
							U_ExtLoc_values[k++] = values[j];
						}
						else
						{
							std::map<int, int>::iterator ret1 = map_neighbor.find(colfrom[j]);
							if (ret1 == map_neighbor.end()) continue;
							std::map<int, int>::iterator ret2 = map_Ext[ret1->second].find(colind[j]);
							if (ret2 == map_Ext[ret1->second].end()) continue;
							U_ExtExt_colind[t] = ret2->second;
							U_ExtExt_values[t++] = values[j];
						}
					}

					delete[] colfrom;
					delete[] colind;
					delete[] values;
				}

				free(recv_buffer[r]);
			}
			else if (recv_ptr[r + 1] > recv_ptr[r])
			{
				for (int i = recv_ptr[r], k = U_ExtLoc_rowptr[i], t = U_ExtExt_rowptr[i]; i < recv_ptr[r + 1]; U_ExtLoc_rowptr[++i] = k, U_ExtExt_rowptr[i] = t)
					;
			}
		}

		delete[] map_Ext;

		int* v = new int[Ext_size];
		double* y = new double[Ext_size];

		for (int i = 0; i < Ext_size; ++i)
			v[i] = -1;

		for (int i = 0; i < n; ++i)
		{
			if (interior[i]) continue;

			w[i] = i;
			x[i] = 0.0;

			for (int j = A_LocLoc_rowptr[i]; j < A_LocLoc_rowptr[i + 1]; ++j)
			{
				int jcol = A_LocLoc_colind[j];
				double jval = A_LocLoc_values[j];

				w[jcol] = i;
				x[jcol] = jval;
			}

			for (int j = A_LocExt_rowptr[i]; j < A_LocExt_rowptr[i + 1]; ++j)
			{
				int jcol = A_LocExt_colind[j];
				double jval = A_LocExt_values[j];

				v[jcol] = i;
				y[jcol] = jval;
			}

			for (int k = L_LocExt_rowptr[i]; k < L_LocExt_rowptr[i + 1]; ++k)
			{
				int kcol = L_LocExt_colind[k];
				double kval = y[kcol] * D_Recip_Ext_values[kcol];

				y[kcol] = kval;

				for (int j = U_ExtExt_rowptr[kcol]; j < U_ExtExt_rowptr[kcol + 1]; ++j)
				{
					int jcol = U_ExtExt_colind[j];
					double jval = U_ExtExt_values[j];

					if (v[jcol] == i)
						y[jcol] -= kval * jval;
				}

				for (int j = U_ExtLoc_rowptr[kcol]; j < U_ExtLoc_rowptr[kcol + 1]; ++j)
				{
					int jcol = U_ExtLoc_colind[j];
					double jval = U_ExtLoc_values[j];

					if (w[jcol] == i)
						x[jcol] -= kval * jval;
				}
			}

			for (int k = L_LocLoc_rowptr[i]; k < L_LocLoc_rowptr[i + 1]; ++k)
			{
				int kcol = L_LocLoc_colind[k];
				double kval = x[kcol] * D_Recip_Loc_values[kcol];

				x[kcol] = kval;

				for (int j = U_LocLoc_rowptr[kcol]; j < U_LocLoc_rowptr[kcol + 1]; ++j)
				{
					int jcol = U_LocLoc_colind[j];
					double jval = U_LocLoc_values[j];

					if (w[jcol] == i)
						x[jcol] -= kval * jval;
				}
			}

			D_Recip_Loc_values[i] = 1.0 / x[i];
			for (int j = L_LocLoc_rowptr[i]; j < L_LocLoc_rowptr[i + 1]; ++j)
				L_LocLoc_values[j] = x[L_LocLoc_colind[j]];
			for (int j = U_LocLoc_rowptr[i]; j < U_LocLoc_rowptr[i + 1]; ++j)
				U_LocLoc_values[j] = x[U_LocLoc_colind[j]];
			for (int j = L_LocExt_rowptr[i]; j < L_LocExt_rowptr[i + 1]; ++j)
				L_LocExt_values[j] = y[L_LocExt_colind[j]];
			for (int j = U_LocExt_rowptr[i]; j < U_LocExt_rowptr[i + 1]; ++j)
				U_LocExt_values[j] = y[U_LocExt_colind[j]];
		}

		delete[] w;
		delete[] x;
		delete[] v;
		delete[] y;

		delete[] D_Recip_Ext_values;
		delete[] U_ExtLoc_rowptr;
		delete[] U_ExtLoc_colind;
		delete[] U_ExtLoc_values;
		delete[] U_ExtExt_rowptr;
		delete[] U_ExtExt_colind;
		delete[] U_ExtExt_values;

		double* D_Recip_Send = A.send_buffer;

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label < neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
			{
				int num_rows = send_ptr[r + 1] - send_ptr[r];
				int num_nonzeros = 0;
				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					num_nonzeros += U_LocLoc_rowptr[send_list[i] + 1] - U_LocLoc_rowptr[send_list[i]] + U_LocExt_rowptr[send_list[i] + 1] - U_LocExt_rowptr[send_list[i]];

				int bufsize = (num_rows + 2 * num_nonzeros) * sizeof(int) + (num_rows + num_nonzeros) * sizeof(double);
				send_buffer[r] = (char*)malloc(bufsize);

				int position = 0;

				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
					D_Recip_Send[i] = D_Recip_Loc_values[send_list[i]];

				MPI_Pack(D_Recip_Send + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_DOUBLE, send_buffer[r], bufsize, &position, comm);

				for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				{
					int irow = send_list[i];
					int length = U_LocLoc_rowptr[irow + 1] - U_LocLoc_rowptr[irow] + U_LocExt_rowptr[irow + 1] - U_LocExt_rowptr[irow];

					int* colfrom = new int[length];
					int* colind = new int[length];
					double* values = new double[length];

					int j = 0;
					for (int k = U_LocLoc_rowptr[irow]; k < U_LocLoc_rowptr[irow + 1]; ++k)
					{
						colfrom[j] = rank;
						colind[j] = U_LocLoc_colind[k];
						values[j++] = U_LocLoc_values[k];
					}
					for (int k = U_LocExt_rowptr[irow]; k < U_LocExt_rowptr[irow + 1]; ++k)
					{
						int kcol = U_LocExt_colind[k];
						colfrom[j] = Ext_ranks[kcol];
						colind[j] = Ext_map[kcol];
						values[j++] = U_LocExt_values[k];
					}

					MPI_Pack(&length, 1, MPI_INT, send_buffer[r], bufsize, &position, comm);
					MPI_Pack(colfrom, length, MPI_INT, send_buffer[r], bufsize, &position, comm);
					MPI_Pack(colind, length, MPI_INT, send_buffer[r], bufsize, &position, comm);
					MPI_Pack(values, length, MPI_DOUBLE, send_buffer[r], bufsize, &position, comm);

					delete[] colfrom;
					delete[] colind;
					delete[] values;
				}

				MPI_Isend(send_buffer[r], bufsize, MPI_PACKED, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
			}
		}

		for (int r = 0; r < num_neighbors; ++r)
		{
			if (label < neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
			{
				MPI_Wait(send_requests + r, &status);
				free(send_buffer[r]);
			}
		}

		delete[] recv_buffer;
		delete[] recv_requests;
		delete[] send_buffer;
		delete[] send_requests;

		delete[] Ext_ranks;
	}
}

int ParCSR_ILU_Preconditioner::InSize() const
{
	return D_Recip_Loc.size;
}

int ParCSR_ILU_Preconditioner::OutSize() const
{
	return D_Recip_Loc.size;
}

void ParCSR_ILU_Preconditioner::Apply(const Par_Vector& b, const Par_Vector& x) const
{
	int n = D_Recip_Loc.size;
	double* D_Recip_Loc_values = D_Recip_Loc.values;
	int* L_LocLoc_rowptr = L_LocLoc.rowptr;
	int* L_LocLoc_colind = L_LocLoc.colind;
	double* L_LocLoc_values = L_LocLoc.values;
	int* L_LocExt_rowptr = L_LocExt.rowptr;
	int* L_LocExt_colind = L_LocExt.colind;
	double* L_LocExt_values = L_LocExt.values;
	int* U_LocLoc_rowptr = U_LocLoc.rowptr;
	int* U_LocLoc_colind = U_LocLoc.colind;
	double* U_LocLoc_values = U_LocLoc.values;
	int* U_LocExt_rowptr = U_LocExt.rowptr;
	int* U_LocExt_colind = U_LocExt.colind;
	double* U_LocExt_values = U_LocExt.values;

	double* bv_Loc = b.Loc.values;
	double* xv_Loc = x.Loc.values;
	double* xv_Ext = x_Ext.values;

	if (n == 0) return;
	
	for (int i = 0; i < n; ++i)
		xv_Loc[i] = bv_Loc[i];

	MPI_Request* recv_requests = new MPI_Request[num_neighbors];
	MPI_Request* send_requests = new MPI_Request[num_neighbors];
	MPI_Status status;

	for (int i = 0; i < n; ++i)
	{
		if (!interior[i]) continue;
		double temp = xv_Loc[i];
		for (int j = L_LocLoc_rowptr[i]; j < L_LocLoc_rowptr[i + 1]; ++j)
			temp -= xv_Loc[L_LocLoc_colind[j]] * L_LocLoc_values[j];
		xv_Loc[i] = temp;
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (label > neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(xv_Ext + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_DOUBLE, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
		if (label > neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int i = 0; i < n; ++i)
	{
		if (interior[i]) continue;
		double temp = xv_Loc[i];
		for (int j = L_LocExt_rowptr[i]; j < L_LocExt_rowptr[i + 1]; ++j)
			temp -= xv_Ext[L_LocExt_colind[j]] * L_LocExt_values[j];
		for (int j = L_LocLoc_rowptr[i]; j < L_LocLoc_rowptr[i + 1]; ++j)
			temp -= xv_Loc[L_LocLoc_colind[j]] * L_LocLoc_values[j];
		xv_Loc[i] = temp;
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (label < neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = xv_Loc[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_DOUBLE, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (label < neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	for (int r = 0; r < num_neighbors; ++r)
		if (label < neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			MPI_Irecv(xv_Ext + recv_ptr[r], recv_ptr[r + 1] - recv_ptr[r], MPI_DOUBLE, neighbor_ranks[r], MPI_TAG, comm, recv_requests + r);

	for (int r = 0; r < num_neighbors; ++r)
		if (label < neighbor_labels[r] && recv_ptr[r + 1] > recv_ptr[r])
			MPI_Wait(recv_requests + r, &status);

	for (int i = n - 1; i >= 0; --i)
	{
		if (interior[i]) continue;
		double temp = xv_Loc[i];
		for (int j = U_LocExt_rowptr[i + 1] - 1; j >= U_LocExt_rowptr[i]; --j)
			temp -= xv_Ext[U_LocExt_colind[j]] * U_LocExt_values[j];
		for (int j = U_LocLoc_rowptr[i + 1] - 1; j >= U_LocLoc_rowptr[i]; --j)
			temp -= xv_Loc[U_LocLoc_colind[j]] * U_LocLoc_values[j];
		xv_Loc[i] = temp * D_Recip_Loc_values[i];
	}

	for (int r = 0; r < num_neighbors; ++r)
	{
		if (label > neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
		{
			for (int i = send_ptr[r]; i < send_ptr[r + 1]; ++i)
				send_buffer[i] = xv_Loc[send_list[i]];
			MPI_Isend(send_buffer + send_ptr[r], send_ptr[r + 1] - send_ptr[r], MPI_DOUBLE, neighbor_ranks[r], MPI_TAG, comm, send_requests + r);
		}
	}

	for (int r = 0; r < num_neighbors; ++r)
		if (label > neighbor_labels[r] && send_ptr[r + 1] > send_ptr[r])
			MPI_Wait(send_requests + r, &status);

	for (int i = n - 1; i >= 0; --i)
	{
		if (!interior[i]) continue;
		double temp = xv_Loc[i];
		for (int j = U_LocLoc_rowptr[i + 1] - 1; j >= U_LocLoc_rowptr[i]; --j)
			temp -= xv_Loc[U_LocLoc_colind[j]] * U_LocLoc_values[j];
		xv_Loc[i] = temp * D_Recip_Loc_values[i];
	}

	delete[] recv_requests;
	delete[] send_requests;
}

}